#include "Fit/RooBreitWignerWave.h"
#include "RooAbsReal.h"
#include "RooFit.h"
#include <iostream>
#include <math.h>

// ClassImp(RooBreitWignerWave);

const double mK      = 0.493677;
const double mpi0    = 0.1349766;
const double mkstar  = 0.88496;
const double mkstar2 = mkstar * mkstar;
const double mphi    = 1.0194550;
const double mphi2   = mphi * mphi;

RooBreitWignerWave::RooBreitWignerWave(const char* name, const char* title, RooAbsReal& _x,
                                       RooAbsReal& _mass, RooAbsReal& _width, RooAbsReal& _l) :
  RooAbsPdf(name, title),
  x("x", "x", this, _x),
  mass("mass", "mass", this, _mass),
  width("width", "width", this, _width),
  l("l", "l", this, _l)
{}

RooBreitWignerWave::RooBreitWignerWave(const RooBreitWignerWave& other, const char* name) :
  RooAbsPdf(other, name),
  x("x", this, other.x),
  mass("mass", this, other.mass),
  width("width", this, other.width),
  l("l", this, other.l)
{}

Double_t RooBreitWignerWave::evaluate() const
{

  double msum2 = (mK + mK);
  double mdif2 = (mK - mK);
  double x2    = x * x;
  msum2 *= msum2;
  mdif2 *= mdif2;

  double a = fabs(x2 - msum2);
  double b = fabs(x2 - mdif2);

  double pK = sqrt(a * b) / (2 * x);

  //************************************obtain c_pK**************
  double tmp_c_pK = (mphi2 - msum2) * (mphi2 - mdif2);
  double c_pK     = sqrt(tmp_c_pK) / (2 * mphi);

  //************************************************************************
  double d_M2 = x * x - mass * mass;
  double coff = pow((pK / c_pK), 2 * l + 1) * (mass / x) * (mass / x);

  double tmp_bw = 1.0 / (d_M2 * d_M2 + coff * coff * mass * mass * width * width);
  double bw     = tmp_bw * mass * width * coff;

  return bw;
}